Mode bifurcation in the Rayleigh- Taylor instability of binary condensates 
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We examine the generation and subsequent evolution of Rayleigh Taylor instability in anisotropic 
binary Bose-Einstein condensates. Considering a pancake-shaped geometry, to initiate the instabil- 
ity we tune the intraspecies interaction and analytically study the normal modes of the interface in 
elliptic cylindrical coordinates. The normal modes are then Mathieu functions and undergoes bifur- 
cation at particular values of anisotropy and ratio of number of atoms. We find that the analytical 
estimates of the bifurcation parameters are in good agreement with the numerical results. 
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I. INTRODUCTION 

Rayleigh- Taylor instability (RTI) [IHS] is the instability 
of an interface between two fluids, which sets in, when 
a layer of lighter fluid supports a denser one or when 
a lighter fluid pushes a denser one, under the influence 
of gravitational field or some external potential. This 
occurs due to unfavourable energy conditions and as a 
result, the fluids tend to swap their positions. Any per- 
turbation arising on the interface, however, small it may 
be, grows exponentially due to RTI and turbulent mix- 
ing of the fluids occur. During the process of mixing, the 
interface gets deformed and develops complicated non- 
linear patterns with mushroom shapes. The phenomenon 
of RTI is widely common in nature, ranging from con- 
vection of water to dusty plasma in atmosphere to su- 
pernova explosions [IHS]. Recently, RTI has also been ob- 
served in a trapped two-species Bose-Einstein condensate 
(TBEC)T], where, intraspecies scattering length plays a 
major role . Systems of trapped TBEC's that have been 
so far studied for observing RTI are a tight, symmetric 
pancake-shaped system in which the components sepa- 
rate out radially, a cigar-shaped trap in which phase- 
separation occurs in the axial direction and a perfectly 
spherical symmetric trap [5141 Oj. Though experimental 
studies on RTI are rare, theoretical studies on interfacial 
instabilities has been a major research topic in the re- 
cent years. Other instabilities such as, Kelvin-Helmholtz 
instability (KHI)[TT], Faraday instability have also been 
predicted in TBEC [HI [T3] . Experimental observation of 
quantum KHI and Faraday waves in BEG can be found 
in Refs.fHIIl]. 

In the present work, we study RTI in a TBEC confined 
in a harmonic trapping potential. The intraspecies and 
interspecies interaction between the atoms are taken to 
repulsive. The initial state of the TBEC that we consider 
for our study, is a phase-separated(immiscible) configu- 
ration in which the species with weaker intraspecies re- 
pulsive interaction is surrounded by the other. In the 



phase-separated domain, the interface of the TBEC is a 
circle when the quasi-two dimensional trap is perfectly 
symmetric. To initiate RTI, we decrease gradually the 
s-wave scattering length of the outer species through a 
magnetic Feshbach resonance. As RTI sets in, the outer 
species tends to sink to the center of the trap and insta- 
bilities begin to occur on the circular interface separating 
the two components. Now, if the anisotropy of the trap 
is increased along a particular direction, the circular in- 
terface evolves into an elliptic cylindrical one. Due to 
RTI, the nature of various non-linear patterns developed 
on the interface changes on varying the geometry of the 
trapping potential. It has been observed that at a critical 
value of the anisotropy parameter, the normal modes on 
the interface bifurcates. 

This paper is organized as follows: In Section [TTj we for- 
mulate the problem using mean-field dynamics in a quasi- 
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two dimensional harmonic trap. In Section [III| i 
discuss about the interface geometry and normal modes 
of the interface and formulate the Helmholtz equation us- 
ing elliptic cylindrical coordinates. In Section |III B[ we 
derive an analytic condition for the temporal decay con- 
stant in an elliptic cylindrical interface. Lastly, in Section 
|IV[ we present numerical results showing the dynamics 
of TBEC as a resuh of RTI. 



II. PHASE SEPARATED PANCAKE SHAPED 
TBECS 

In the mean field approximation, the TBEC is de- 
scribed by a set of coupled Gross-Pitaevskii equations 
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where i — 1, 2 is the species index, ga — Anh^ai/mi with 
rrii as mass and a,; as s-wave scattering length, is the 
intra-species interaction; gij = 2'Kh? aij / niij with rriij = 
mirrij /{rrii+mj) as reduced mass and Uij as inter-species 
scattering length, is inter-species interaction and ^i is 
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the chemical potential of the zth species. The trapping 
potential is 



(2) 



where, uj is the radial trap frequency, considered identi- 
cal, for the two components, and ai , Xi are the anisotropy 
parameters. For simplicity of analysis, we consider trap 
potentials of both the species have the same geometry 
tti = Q;2 = a , Ai = A2 = A and mi — TO2 = m. The 
energy of the TBEC is 
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To express the energy in suitable units, we define the os- 
cillator length of the trapping potential Oosc = i/^/ (ww) 
and consider huj as the unit of energy. We then divide 
Eq.([3| by hiu and apply the transformations x = x/aosc 
, y = y/oosc, z = z/flosc, i = tuj, and E = E/{huj). The 
transformed order parameter 
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Using Thomas-Fermi approximation in 
Eq.(|8]), one can show the two components are phase- 
separated when {M\2 > \/A/iiA/'22), where, Afu and A/ij 
are all positive. To examine RTI, we consider the phase 
separated state in axis symmetric trapping potentials 
with coincident centers and numerically solve the pair 
of time-dependent GP equations 



vl + v,{x,y) + J2^f^A^l{^,yW 



to study the dynamical evolution. 
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and energy of TBEC in scaled units is given by 
2 
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where, Uu = inau/aosc and U12 = 4T:ai2/aosc- For 
simplicity of notations, from here on we will represent 
the transformed quantities without tilde. Thus, in scaled 
units, the coupled 3D GP equation is given by 



V' + V,ix,y,z)+J2G^A^A 



(6) 



where, Gu — NiUu and Gij — NjUij. For the present 
work, we consider a pancake shaped trap, the axial fre- 
quency is much larger than the radial frequency(A ^ 
1). In this situation, the transformed order parameter 
$(a;, y, z) is factorized into 



$(x,y,z) = 0(x,y)C(2:) 
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where, ^(z) is the normalized state of axial trapping po- 
tential V,"""^"^ = A^z^/2. From Eq.([6| after integrating 
out the axial order parameter, we obtain the scaled cou- 
pled 2D GP equations 



-Vi + V:,(a;,y)+^M,|0,(x,2/)p 



(8) 



FIG. 1. Phase separated profiles of ^^Rb-^'^Rb mixture at 
t = 0. The figure on the left shows the inner species(*^Rb) 
with 022 ~ 99aB. The figure on the right shows the outer 
species(*^Rb) with an = 460aB. 



III. INTERFACE GEOMETRY AND MODES 

In the phase-separated domain, the interface of the 
TBEC is a circle when a is unity. It is, however, trans- 
formed to an ellipse when a > 1 . A typical density profile 
of the the phase-separated TBEC with a = 1 is shown in 
Fig. [1] Compared to Eq. (pj), a more general form of 2D 
trapping potential is V{x,y)— m{uj'^x'^ +ujyy'^)/2, where 



represent angular trapping frequency along x or y. 



Defining the geometric mean oj 
potential is 

1 



y/uJxU)y, the trapping 



y 



r 



(10) 



x 

where, /3 — Cu/ojx = \J a and 7 
The density distribution of the TBEC, at moderate 
anisotropics, follows the geometry of the trapping poten- 
tial. At larger anisotropics the interface energies modi- 
fies the density distribution and leads to difference from 
the geometry of the trapping potential. For the present 
study, we consider the TBEC at moderate anisotropies. 
The interface of the TBEC is then an ellipse 
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corresponding to the anisotropy parameters of the trap- 
ping potential. 
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FIG. 2. Phase separated profiles of *^Rb-*^Rb mixture at 
t = 0. The figure on tfie left shows the inner species(*^Rb) 
with a22 ~ 99aB. The figure on the right shows the outer 
species(**Rb) with an = 460aB. 

At the interface, the densities are low, neglecting the 
intraspecies and interspecies interaction term, we get 
from Eq.([8| 

(-vi + v,)(i)^ = /I, 

2 

where, /I^ = ^, - ^Afij\(j)j{x,y)\'^ . 
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Using Eq.(ll) and Eq.([l2|), we get (in scaled units) 

vi(/),+ ( /I. = 0- (13) 



Defining the parameter fc| = /I^ — a, the equation is 



= 0. 



(14) 



This is the Helmholtz equation in 2D. It must, however, 
be emphasized that the equation is valid only at the inter- 
face or close to it. Away from the interface the densities 
are not small and intraspecies interactions is large. 



Using separation of variables, the equation is simplified 
to the Mathieu equations [HI [T7] 
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where, A ^ A + and q = a^k'^/A. Here A is the 

separation constant and returning to the earlier definition 
of the trapping potential, the anisotropy parameter a — 
(3/ J. The interface is an ellipse with eccentricity e — 
\/l ^ 1/a^ and from the theory of conic sections a ~ 
j3e — \foLe. Based on these definitions, the constants in 



the Eq. (17) and (18) are redefined as 
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The constants in this form are easier to connect with the 



parameters of trapping potentials. The Eqns.( 17 and 18 ) 
then assumes the form 
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The interface of the TBEC, an ellipse, has fixed coordi- 
nate u representing the elliptic cylinder. But the angle 
coordinate v varies and lies in the domain [0, 27r). Thus 
0, solutions of the second equation, represent the normal 
modes of the interface. For circular interface, a = 1 and 
e = 0, only < A is physically admissible and the solu- 
tion of the equation is reduced to sinusoidal functions. 



Normal modes of the interface 



B. Instability at the interface 



For linear stability analysis of the interface modes due 
to a small perturbation, to identify the onset of RTI, we 
transform the Eq. ( 14 1 to elliptic cylindrical coordinates 
{u,v). Here, the coordinate v represent the asymptotic 
angle of confocal hyperbolic cylinders symmetrical about 
the X-axis. And, the u coordinates are confocal elliptic 
cylinders centered on the origin[16j. The transformation 
is defined by the relations x = a cosh m cost), and y = 
asinhitsinu and a is the focal distance along cc-axis. We 
take the coordinates on the z = plane as we consider 
the TBEC in 2D. The Eq. (IT4| then assumes the form 
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where (/) is the solution of the form </> — U{u)Q{v). Sub- 
stituting, (j) back in Eq.(15l we get, 

/I d^U 



U dv? 
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e dv^ 



+ c^sin^w = 0. (16) 



For TBEC in traps, the gradient of the trapping po- 
tential is like the gravitational force in the conventional 
fluid dynamics and the flows within TBEC is modelled 
as potential flows. Consider the interface of the TBEC, 
using the method of normal modes, any arbitrary distur- 
bance on the interface may be resolved into independent 
modes of the form 



(23) 
(24) 
(25) 



Here, ^ is the position of the interface relative to the 
equilibrium configuration, and (p^ is the increments in the 
velocity potential of the ith species about the interfacial 
region caused due to disturbance in the system. ^ and 
(pi are the amplitude of the modes and s is the temporal 
decay constant. 
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We know that for an incompressible fluid, the Lapla- 
cian of the velocity potential vanishes V^(/)j = and from 
the expression of the normal modes 
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Using separation of variables, Eq.(26) can be simplified 
to. 



0. 



(27) 



where, = -{\IQ){d'^Q j dv^). The general solution of 
the above equation in the regions of the two species are 



</>i(u) = ^le 
?2e 



-Cxi 



I2 



Cu 



(28) 



where, Ai, Bi, A2, B2 are arbitrary constants. Here, 
it is to be noted that the sign of the exponents in the 
two solutions are interchanged. This is to indicate that 
the relative distance from the interface, within the two 
species, are in opposite directions. We may recall that 
the instabilities occur only at the interface or close to it. 
At any point far removed from the interface the normal 
modes must decay to zero. Thus, normal modes are of 
the form (jfi = Ai exp(— Cu) and 02 = ^2 exp(CM). The 
velocity potentials in Eq. ( 24 ) and Eq. ( 25 ) are then 



(29) 
(30) 



The dynamical evolution of the interface is described 
through a combination of the continuity equation, Eu- 
ler's equation and Bernoulli's theorem [3l I18j. For stabil- 
ity analysis of the interface, we linearize these equations 
and neglect quadratic terms in 0]^, 02 ^-iid ^. After lin- 
earization, we get 
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where, g is the gradient of the trapping potential V{x, y). 
On the interface, using the solutions obtained earlier, 
from Eq. ( 32 ) one can show that 



ni{sAi + g£,) = n2 (s^2 + gS.) ■ (33) 
and A2 ~ s^/C. Using these values in the above equation 



In a similar way, from Eq. (311, we obtain Ai 



ni 



C 
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Simplifying this equation, one arrives at the definition of 
the temporal decay constant 



± 



Cg{ni - na) 
ni + n2 



(35) 



The densities of the condensates ni and 712 are at a point 
(m, v) on the interface. We recollect that 712 refer to the 
density of species at the center which is surrounded by 
the species with density ni . The interface is stable when 
s is imaginary ( ni < 712) and oscillates when perturbed. 
However, when ni > 7^2, the value of s is real and any 
perturbation, however small, grows exponentially with 
time. This is the prerequisite for RTI in binary conden- 
sates. In this context, Atwood number F is given by [TH] 



iFrom Eq. (22 I 



ni - 712 

Til -t- 7l2 ' 



A + ^^k^e^a (1 - cos2w) 



(36) 



(37) 



the temporal decay constant is then 

1 

s = ±[A + 2g(l-cos2i;)]J ^ (38) 

[ 711+77-2 J 

where for compact notation we ha ve u sed the relation 



given earlier. Thus Eq. 

as 



38 



can be rewritten 



s = ± + 2g (1 - cos 2-1;)] * yJVg. 



(39) 



This shows that s is a function of v, the angular coordi- 
nate. 




FIG. 3. Shaded regions indicate the values in aq-plane where 
solutions of angular Mathieu equation, ce^iv, q) and se^(7;, q), 
exist. The black (red) colored curves are the values of a and 
q for which integer order, ce„_i(u,g) (se„(t;,g) ) with n = 
1,2,3, . . ., solutions exist. 
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Allowed solutions 



The solutions of Eq. ( 22 ) , the angular Mathieu equa- 
tion, are the cei,(i', g) and se^{v,q) functions |20], cosine 
and sine elliptic functions, respectively. Here, i' is real 
number and denotes the order of the elliptic functions. 
The solutions, however, exist only for certain range of 
a and q, and these are shown as shaded regions in Fig. 
|3j In the figure, the shaded region consists of lobes and 
each are bounded by elliptic function of integer orders 
ce„_i(w,g) and se„(u,g), where n = 1,2,3,.... For the 
present case, when RTI sets in, the mushroom shaped 
superfluid flows have four fold symmetry in the case 
of circular symmetry. So that the flow retains symme- 
try or shape invariance along perpendicular directions. 
The corresponding solution of Eq. (22) which satisfy 



this condition is then ce2 («,(?), and it has the proper- 
ties ce2{v,q) = ce2(v + Tr,q) and ce2{v,q) = ce2{—v,q). 
The loci of the a-q pairings which allow this solution is 
the labeled curve in Fig. [3] 

One property of ce2{v,q) is, the maximum at f = 
undergoes a smooth bifurcation at higher values of q. 
Coming to the description of the interface in the bi- 
nary condensates, from Eq. (19), q is a linear function 
of the anisotropy parameter a. So, as we increase the 
anisotropy the mushroom shaped flows in RTI must un- 
dergo bifurcation. At some value of a, instead of four 
there must be six mushroom shaped inward superfluid 
flow. 



over real-time, we gradually change the scattering length 
over time steps, and study its dynamics [8j. Density pro- 
files as shown in Fig. [4] are the numerical solutions ob- 
tained by this method. 
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FIG. 4. (a)-(c) shows development of mushroom shape pat- 
tern on the interface after t = 358 ms, t = 378 ms, t — 400 
ms. The scattering length is decreased from an = 460ob to 
ail = 55aB between t — ms and t = 200 ms , after that an 
is fixed to 55aB upto t — 400 ms. The images on the upper 
panel correspond to the inner species (*^Rb) and the images 
on the lower panel correspond to the outer species (**^Rb). 



A. Mode bifurcation and density profiles 



IV. NUMERICAL RESULTS 

To corroborate the analytic results for the interface 
modes, as mentioned earlier, we numerically solve the 
pair of coupled Eq^^ We resort to split-step Crank- 
Nicholson method[2T] implemented for binary conden- 
sates. We discretize Eq. [9] both in space and time, and 
propagate the resulting discretized equation in imaginary 
time, over small time steps. In imaginary-time propa- 
gation method, t in Eq. [9] is replaced by —it. This 
method seems to be more appropriate as the stationary 
ground state wave function of the TBEC is essentially 
real and dealing with real variables is more convenient 
than imaginary ones. The split-step imaginary time so- 
lution obtained in a self-consistent way after several it- 
erations, is the stationary state of TBEC for the given 
parameters used in this paper. The time-independent 
solution of TBEC, thus obtained, is used as an initial 
state for real-time propagation. The real-time propaga- 
tion method yields the solution of time-dependent CP 
equation for TBEC, which is used to study the dynami- 
cal behaviour of TBEC. 

As a representative case, we numerically calculate the 
stationary state solution of TBEC based on the afore- 
mentioned method, with the parameters given in Ref. 
[22j . To study RTI, we use the imaginary-time solution 
as the initial state. With the propagation of this solution 



We consider a system of Rb- Rb atoms in a sym- 
metric 2-D harmonic trapping potential with = 
27r X (8, 90)Hz. We choose initial state to be the ground 
state for which an = 460aB, 022 = 99aB, ai2 = 121 = 
214aB, with ob being the Bohr radius. The number of 
atoms are iVi = 5 x 10^ and N2 = 10^ [22 . In this config- 
uration, component 1 (outer), ^'"'Rb completely surrounds 
component 2 (inner), *^Rb. Fig. [I] shows phase separated 
profiles of the TBEC at t = in a perfectly symmetric 
pancake shaped trap i.e. a = 1. 

Now, the s-wave scattering length an of the outer 
species is decreased gradually over time, experimentally 
this is possible through the ^^Rb-'^^Rb magnetic Fesh- 
bach resonance}^. However, throughout the process, 
we maintain (A/12 > V ■^11-^22) so that the TBEC re- 
mains in the immiscible domain. A stage is reached 
when an < 022, ^^Rb-®^Rb interaction weaker than the 
^^Rb-^^Rb interaction. In this situation, the existing spa- 
tial structure of the system is energetically unfavourable 
and the outer species starts penetrating inside the in- 
ner species. Instabilities begin to occur at the interface 
of the two components and eventually grows into a four 
fold mushroom shape superfluid flow as shown in Fig. [4] 

The dynamics and the formation of lobes also depends 
on the geometry of the interface. As the anisotropy of 
the trap a is increased keeping A and other remaining 
parameters fixed, the circular interface evolves into an 
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FIG. 5. (a)-(d) Development of various non linear patterns 
for a = 1, a = 1.2, a — 1.4 and a = 1.6. The images on 
the upper panel correspond to the inner species (*^Rb) and 
the images on the lower panel correspond to the outer species 
(^^^Rb) 



elliptic interface. The penetration of the heavier fluid 
into the lighter fluid gets initiated along the x-axis, fol- 
lowed by the formation of lobes. This happens because, 
the interface is more curved along this direction with less 
confinement. Larger is the curvature, higher is the rate 
of inflow of the heavier fluid. Mass transport gradually 
occurs along y-direction, which is tightly conflned. The 
interface here, is relatively flat and the lobes are formed 
at later stages of evolution. This is clearly evident from 
the superfluid flow pattern soon after the onset of RTI as 
shown in Fig. [5] In Fig. [5j;-d, the lobes of 712 along the 
X-axis are well developed and located deep within ni. As 
a is increased further, the interface is deformed further. 
The lobe along the y-axis undergoes a bifurcation when 
the anisotropy is such that a = 3 and the density proflle 
of the superfluid flow is shown in Fig. [6] Taking an aver- 
age along the interface and close to the bulk of ni, when 
the mode bifurcates the value of /i is 5.25. This can be 
related qualitatively to the analytic results, in which case 
ce2{v,q) undergoes bifurcation at around q « 3.7. Thus 
our numerical results is in agreement with the inferences 
drawn from the analytic solutions of the interface modes. 
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FIG. 6. (a)-(c) Development of various non linear patterns 
for a = 1.8, Q = 2.0 &i a — 3.0. The images on the upper 
panel correspond to the inner species (*^Rb) and the images 
on the lower panel correspond to the outer species (*^Rb) 



B. Bogoliubov analysis 

For a more detailed understanding of the instabil- 
ity, we perform a Bogoliubov analysis for TBEC in a 
2-D harmonic trap. Setting, $i = 4>i + 6(j)i, we ex- 
pand the set of coupled GP equations in Eq. |9] in 
S(j)i{x, y). Here, 0i's fixes the condensate density through 
ni{x,y) = \(j)i{x,y)\^ and 5(j)iS are the deviations from 
the initial ground state, which includes the quasi-particle 
excitations. We consider excitation mode of the form 



(40) 



where, /i^ is the chemical potential, w is the excitation 
frequency, and Ui and Vi are the Bogoliubov amplitudes. 
Using this ansatz, the Bogoliubov equations are 



2m; 
2m. 



■Vi + Vi + 2n.iUii + UijUj - 



+ Uij^/r^nj{uj -Vj) = hwui, (41) 
V]_ + Vi + 2niUii + UijUj ~ fii Vi - UiiTiiUi 

+ Uijy/ninj[vj — Uj) — —hujVi.{42) 



These equations are then numerically diagonalized to cal- 
culate the excitation spectrum. If the frequencies are 
real, the perturbations remain bounded and the system 
is dynamically stable. On the other hand, pure imag- 
inary eigenfrequencies denote instability of the system. 
The eigenmode corresponding to this complex frequency 
grows exponentially and is a signature of dynamically 
unstable system [24] . 

As a case study, we choose N2/N1 = 2. For the 
isotropic case, a = 1, we expand the Bogoliubov ampli- 
tudes in harmonic oscillator basis wave function. When 
(5(/>i is small, Re(a;) increases monotonically upto a critical 
point as the an is decreased. When an < 022 and RTI 
sets in, the low lying excitation modes w starts becom- 
ing imaginary. The value of the Im(w) increases mono- 
tonically as an is decreased further and away from the 
critical point. These imaginary modes are signatures of 
instability in the dynamics of the binary condensates. 



C. Effect of noise 

Numerical studies that have been carried out so far 
are at zero temperature and without any imperfections, 
hence quite ideal. But in experiments, conditions are far 
from ideal. Fluctuations play a major role, and if large, 
may destroy the observed signatures predicted from the 
numerical simulations. One immediate remedy is to in- 
clude fluctuation to our calculations. We introduce white 
noise during the real time evolution of TBEC. The white 
noise is at the level of 0.01%. Even after introducing 
noise, we still observe signatures of RTI as a result of 
changing an. The thermodynamical quantities such as 
energy, chemical potential may vary quantitatively, but. 
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there is no qualitative difference in tlie shape of the inter- 
face after RTI is initiated. Bifurcation of normal modes 
on the interface are still observed at the predicted values 
of the anisotropy of the trap. 

V. CONCLUSIONS 

We have examined RTI at the interface of binary con- 
densates as a function of anisotropy parameter and ra- 
tio of number of atoms. The mushroom shaped super- 
fluid flow is four lobed, as expected, at low anisotropies. 
Based on the analytical studies, the lowest natural mode 
is cen-i{v, q), which describes the four lobed superfluid 



flow. However, at higher anisotropies corresponding to 
larger values of q, one of the maxima of ce2(v,q) bifur- 
cates. This is also observed in the numerical simulation 
of the RTI at higher a. The RTI and bifurcation of the 
mode is robust, and observable in presence of white noise. 
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